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Ai3traot - This paper presents a recursive algorithmic 
implementation of the prewindowed high performance 
method of ARMA spectral modeling as described in Part 1. 
This algorithm provides updates of the ARMA models 
optimal autoregressive coefficients (in actuality pre¬ 
diction errors) as each new data point becomes 
available. The algorithm is computationally fast in 
the sense that it requires 0(p) multiplications and 
additions for each update. It is shown that this fast 
recursive algorithm may be implemented using a lattice 
filter arrangement, and it therefore exhibits several 
of the "nice” properties associated with lattice type 
algorithms such as numerical robustness and good con- 
re rgence properties. v 

I. INTRODUCTION 

In Part 1 of this paper we described an algorithm for 
ibtalning an estimate of the power spectral density 
.ssociated with a given time series (x(n) f. In parti- 
ular. the following ARMA (p,q) spectral density model 
as hypothesized 
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A closed form algorithm for estimating this model's 
a* and bj parameters was then developed in which the 
finite set of time series observations 


*(1>, x(2),. 


x(n) 


( 2 ) 


were used in the parameter selection process. This so 
called "high performance” algorithm operates on a data 
block of length n to obtain the model's coefficients 
in a single computational effort. It is therefore 
called a block processing algorithm. 

There are many situations, however, in which a block 
processing algorithm for spectral estimation is not an 
appropriate tool. In a variety of applications, the data 
measurement is an ongoing process and it is therefore 
desirable to recursively update the autoregressive and 
moving average parameters as each new data point 
becomes available. This capability is of particular 
importance in chose cases where one wishes to adaptive¬ 
ly model the spectrum of a long, ongoing time series. 
Algorithms with this recursive updating capability are 
called "recursive algorithms”. 

Recently, several fast recursive spectral estimation 
algorithms have been developed [11—[7]. Most of these 
algorithms are based on so called least-squares AR 
spectral estimation methods [81, although a few ARMA 
algorithms have also been developed [4] i [6], These 
algorithms saek to minimize a prediction error vector 
in order to obtain the desired spectral estimates. 

In this paper we develop a more effective ARMA 
algorithm that efficiently computes Che optimal auto- 
ragraaslve coefficients by recursively updating a sec 
of prediction error elements as each new data point la 
obaerved. This recursive algorithm is based on the 
prewindowed high performance method as described in 
Part 1, and, therefore is predicated on che approxi¬ 
mation of the ARMA model's underlying Yule-Walker 
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equations. Although similar algorithms may be derived 
for che unmodified and for che two other modified 
versions of che high performance method, the recursive 
algorithm based on che prewindow version is a blc 
easier to derive and is characterized by a "fast start¬ 
up" capability in that spectral estimates are poaslble 
with as few as two data points. The recursive algorithm 
herein presented is computationally efficient in the 
sense chat 0(p) multiplications and additions are 
required to update the "necessary" parameters as each 
new data point is observed. This paper's recursive 
algorithm was originally developed in [9). A more 
straightforward derivation is herein presented which 
provides a greater degree of insight. In addition, a 
lattice filter implementation of this algorithm is 
developed. Moreover, because of this ladder-type 
Implementation, chis algorithm is characterized by 
several other nice properties associated with ladder 
algorithms such as numerical robustness and good 
convergence properties [10], [11]. 

II. THE PREDICTION ERROR VECTORS 

The recursive update equations herein presented do 
not explicitly update che ARMA model's autoregressive 
coefficients in obtaining optimal updated spectral 
estimates. Instead, a set of "equivalent" parameters 
known as prediction errors are updated. In this 
section, we discuss the relationship between the pre¬ 
diction errors and autoregressive coefficients. 

As outlined in Part 1 of this paper, the optimal p ch 
order set of autoregressive coefficients for the pre¬ 
windowed version of the high performance method are 
obtained by solving the following system of p linear 
equations in p unknowns 


[Y + X la + Y + X • 8 

n,p n,p —p n,p —ti — 

(3) 

where 


x • (x(l), x(2),. . x(a>] 

(4a) 


(4b) 
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in which 8 is the zero vector and S is the down shift 
operator.T Here the dagger symbol (+) denotes complex 
conjugate transposition and the prime symbol (') denotes 

transposition. The subscripts p and n explicitly , - 

indicate that the denominator order of the spectral 
model of equation (1) is p and thac n data points are 
available. Whenever this explicit information is not 
needed, we will use x, £, X, and Y in place of Xj,, ^ 1> , 

X n p , and Y n p, respectively. ’ 

it is recalled from Part 1, that the high performance 


ARMA modeling approach is predicated on approximating t 

''In this paper, matrices are denoted by capital letters 
(e.g. X), vectors are denoted by underlined lower case 
English letters (e.g. x) and scalars are denoted by 
lower case Greek letters (e.g. a). Moreover, the down 
shift operator S Is defined by 
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Yule-Walker equation* where t>p. Upon examination of 
expression (3), it is apparent chat we have hare re¬ 
stricted t"p. This restriction is required in order 
to facilitate cha development of the fast recursive 
algorithm. Unfortunately, by requiring t • p, the 
spectral eatlaation performance suffers in comparison 
to chat achieved with larger values of c. as the data 
length n increases, however, this performance degrad¬ 
ation diminishes and typically is of an insignificant 
nature. This is indeed fortunate since it is precisely 
for long data length cases that cha recursive algorithm 
would most likely be utilized. 

Under the assumption chat Y+X is nonsingular, Che 
optimal autoregressive coefficient vector which satis¬ 
fies expression (3) is given by 

a • - -[Yt X„ ]'V x (5) 

—p 1 n,p n,p' n,p-n 

In what is co follow, it is beneficial Co Interpret 
this autocorrelation coefficient selection procedure 
from a prediction error viewpoint. Namely, we may re¬ 
formulate expression (3) as 

Y + f x - 8 (6) 

n.p —p,n - 

in which £p X n is the so-called forward prediction error 
as specified by 


It is referred to as the forward prediction error since 
Its kth component can be interpreted as being the error 
resulting from a prediction of the element x(k) by a 
linear combination of the p most recent time series 
elements x(k-l), x(k-2),. . ., x(k-p). 

The optimal autoregressive coefficient vector (5) can 
be then associated with an auxiliary minimization 
problem involving the prediction error vector. Namely, 
ic is readily shown chat this optimal vector minimizes 
Che following quadratic functional 

g(a ) - f* *W f* n (8) 

—p - p, n —p, n 

where W Is Che n*n positive semideflnlte matrix 
specified by 

W - Y Y * (9) 

n.p n.p 

To reinforce this prediction error interpretstion, let 
us define che following estimate of vector Xg 

x - -X a (10) 

n n.p —p 

which in turn generates the forward predicting error 

JL*„ * * ' L (U> 

"■p«n ~o 

Upon substitution of expression (5) into (10) the opti¬ 
mum forward prediction error vector is given by 


in’ ' Vp^ntpVp 1 " 1 ^ 2n 


while che minimizing forward prediction error for this 
selection becomes 

Ip.n- (* - X n.pt Y ntp X n.pi'\‘pK < 12 » 

■ P c x 
XY in 

We have here used the compact matrix product represent¬ 
ations 4 _ 

P„ - X„ [Y * X ]Y n (13a) 

XY n,p' n.p n.p n.p 


P c ■ I - P 
XY XY 

Since we are only interested In th* optimal Xg and 

f x , w* will drop th# symbol and asauM that x 
» D “ 

and f x ar# tha optimal ona* a* given by equation (12) 

P»n 

W« may alto daftna cha dalayad backward prediction 
error vaccor for *g by 2 

d x - S^x + X n i (14) 

■T>,n -n n,p —p 

It can be seen that the k c>> row of aquation (14) 
represents a prediction of x(k-p-l) by a linear combina¬ 
tion of the p moat immediate future valuas x(k-p), 
x(k-p+l).x(k-l). Th* resulting error in this 

backward prediction is d x (k). In this case the 
P> n 

optimum *p vector is the ona that minimizes th* 
quadratic function 

8(i ) * [d* n l t Wd* (15) 

—P “p,n p.n 

where w Is defined in equation (9) 

In this case che optimal a la given by 

a* - -(Y f X ] _1 Y + (S 1 ^ 1 * ) (16) 

P n,p n,p n.p —n 

In a similar manner to the forward prediction error case 
the optimal estimate of sP^Xg is specified by 


-X a 
n,P 


’ p XY< sP+1 Sn> 


and the optimal delayed backward prediction error vector 

x .. s p+l x _ S^x*- p c (S^x ) (17b) 

—p . n -Ti —n a i —n 

where Pjy and Pyy are given by equation^(13). Further¬ 
more, It can be shown that the optimal a p as given by 
equation (16) also arises by approximating p Yule- 
Ualker equations In a manner similar to the approxima¬ 
tion given by equation (3) for 

It is clear that the forward prediction error vector 

fx and cha autoregressive coefficient vector *_ are 
■^p.n —P 

interchangeable in che sense chat one can always be 
found from che ocher using equation (7). Similarly 
d x and ip are interchangeable since ona can always 
b8’?ound from tbs ocher using equadon (14). It is 
also true that the 2p elements of Sp and ip are Inter¬ 
changeable with the 2p elements f x Q (n), fj n (n),..., 

f x (n) end d* (n), d* (n), ..., d x (n) (i.*., cha 
P.n l.n z.n p.n 

nth elements of the 2p prediction error vactors £ x g , 

f? .f x and d? , d? .d x ). 

—2,n ' -p, n —l.n —2,n P.n 

We will show this last fact In Section VIII, where 
we will also see that the prediction errors lead to a 
lattice filter structure which is related co the auto¬ 
regressive coefficient vectors. 

In the fast recursive algorithm, the autoregressive 
coefficient vectors £p and ip are not directly updated. 

Instead, the prediction error clamant* f x (n). 

i.n 

f x (n), and d x (n),..., d x (n) are updated. Sine* 
p.n i.n p.n 

these elements are Interchangeable with the auto¬ 
regressive coefficients, there is no information lost 
in updating only the prediction error elements. How¬ 
ever, che prediction error elements may ba updated In 
a computationally efficient manner, requiring O(p) 
multiplication* and additions for th* update. More¬ 
over, at we shall see later, Che 2p prediction error 
elements enable us to find all of tha Sg and *„, 
vectors for ARMA denominator orders from 1 to p. It 
Is for these reasons that we ch oose to update th* 

^We use the tern "delayed" because although th# subscript 
n appears, x(n) is never used in (14). The undelayed back¬ 
ward prediction error vector will be discussed in a later 
section. 
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prediction error elements. 

He may also obtain additional autoregreeaive co¬ 
efficient estimates similar to a and a by considering 
the prediction error vectors associated with the vector 
given In equation (4b). Specifically, the opclaaa 
forward prediction error corresponding to is defined 
as 

f p • v -*• Y c* (18a) 

-p,n n,p =p 

where 

C* ' -f*n + n In <18b) 

—p ii|p a,p Q,p ni 

end X n ,p tn.p vui fn ar * defined in equation (4). As 
in earlier cases, c£ can be found by approximating p 
Yule-Walker equations or, equivalently, by minimizing 
the quadratic functional 

h <V * <0 +[1 n,p ^IpHl^.J < 19 > 

Corresponding to this optimal autoregraaalve co¬ 
efficient vector we may also define the estimate 

in*' - Y n. P ^ ■ ^.p^ntp^.p 1 ' 1 ^!^ 

- P Y X ^ < 20 > 

which, in turn, gives rise to che optimum error vector 


where P^ x and P^ are defined as in equation (13) 

Finally, the opclmal delayed backward prediction 
error vector for ^ is defined by 


d y S^v + Y c 

Sp.n *n ti,p -p 


V • - tV.p 1t n.pl'Vp (Spfl Ai > 


Again, w* can define the predicted value of Sjr by 

S ^\‘* 1 f n.pt*nVn f pl"Vp <Spfl *> * f YX (SP+1 i ) 

( 22 ) 

and it follows that 

4.„‘- ,st *‘v m> 

Just as for the Xn vector, the 2 p entities 

.*£,„<•>• and ^.tv (n> '--’ ^,n (n) C " be 

efficiently updated and enable ua to determine the 
optimal c,,, and Cj, coefficient vectors for all ARMA 
model denominator orders from 1 to p.2 


III. THE HILBERT SPACE SETTING 

The problem of recursively updating tha prediction 
error vectors in the fast algorithm can bs more easily 
understood by casting the problem in a Hilbert space 
setting. Consider che n dimensional complex 
Euclidian space 


R - C" 


1 c *... *c 


with che standard vector inner produce defined by 
n 

t r * 

<x,i> * j i * I *(1) y(l) 

1-1 


’Throughout the remainder of the paper, the symbol 
will be dropped and the prediction error vectors 
fX, J x , fy, and dY are assumed to be the optimal ones. 


He note that the nxl vectors Xg, S**„, y n , and S*^ axe 
all elements of 8. Moreover, the p columns of matrix 
Xq. arc also alamenta of B. Tha mat of all linear 
combinations of these p e le m en ts la a subepace of B. 
which wa denote by Mg. Similarly, My is the subepace 
■panned by tb* p columna of Y n< p. 

Let ua now consider tha forward prediction of x^. From 
aquation (12a) we tea that Xgla formed by a matrix 
multiplication involving xp. The matrix Pgv la seen to 
be a linear operator on tha Hilbert space a. It is ap¬ 
parent that ?jy maps alamenta of B into alamenta la tha 
subapace Mg, that la 

P *Y= * ~ «X . (26) 

Also, It is evident from aquation (13) that Pgy • Pjy 
so chat the operator Pgy la a projection operator onto 
the subapace Mg. In general, Pgy is not tha orthogonal 
projection operator onto subapace Mg. Instead, tha as¬ 
sociated direction of projection la determined by che 
matrix Y n .. It can ba seen from aquation (6) that tha 
dlrectlon’of projection of Pgy la orthogonal to My. 
Thus, Pgy is the projection operator onto che subspace 
Mg along My4 (the orthogonal complement of My). 

With those thoughts in mind, wa can provide a simple 
geometric Interpretation to tha four error vectors 
described in the last secclon. In particular, the 
geometric relationship between Xj,, x„, and f* Q la 
depicted in Figure 1. ’ 



Figure 1: Geometric Relationship Between x n , x 0 , “8 
the optimal prediction error is £*~ n 

The vector is seen to be thet projection of j^onto 
Kg that Is orthogonal to My. He note from Figure 1 chat 

Jp,ni«y < 27 *> 

or, equivalently, chat 


-P.o’ S 2-n > 


1,2 .. 


The geometric relationships for d* a , fX n , and d£ n *re 
similar to Figure 1. 

Since Pjy and Pyg are projection operators, so sra 
their complements P^y and P$g. It follows that 

(?L ) 2 • PL (28a) 


He can also saa from aquation (13) that P XY and P YX 
ars strongly rslacsd, namely 

P XY * t P Y X J + (29 *> 

and P^ - IP^ X )* (29b) 

In addition to tha four prsdictlon error vectors, 
there ars four inner products thet are useful in deriving 
the fast algorithm. Theae complex-valued scalert ere 
defined as: 

°p,n 4 ‘ipVXn 1 ’ ^ )tp n lS ^2n ! < 30) 
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r i (d X ] T [f y 1 
p,n — 1 p t n* L “P»Q J 


(S^xJ 


,n ^ " ISa^nto.! 


t„c 


v P ,n 4 ^.n^.n 1 ' 

IV. THE PROJECTION OPERATOR THEOREM 


cf** l+» c f«P +1 „ 


(31) 

(32) 

(33) 


From Che results of che laic section lc is apparent 
chat che various prediction error veccors and scalars 
are all described by the operators Pgy and Pyx- As a 
new data point x(n+l) becomes available, we desire to 
update the prediction errors vectors and scalars In a 
computationally efficient manner. Because che opera¬ 
tors ?xy and Pyx are use< 3 repeatedly and their struct¬ 
ures change as new data points become available, we 
prefer co update them and then obtain updated error 
vectors by applying these updated projection operators. 
Recursive update equations for the projection operators 
Pyy and Pyx, are readily obtained by appealing to the 
following theorem. 

Theorem 1. (Projection Operacor Theorem) Let A and 
B be n*m matrices. Furthermore, consider Che aug¬ 
mented matrices A - [A ( a] and B « (B I bj in which 
a and are n*l vectors. If [A *B]“^ and [A + Bl - ^ 
exist, chen the associated projection operator corres¬ 


ponding co the augmented matrices Is given by 


P X5 

* 4 - P AB it'V P L 

(34) 

where 

P AB - AtB'Aj^B^ 

(35a) 


4 * 1 - P AB 

(35b) 


_The theorem may be straightforwardly proven by writing 
[A “B] In cerms of Schur complements and performing 
some matrix algebra. Alternatively, che theorem may 
be proven using Hilbert space concepts. 

L'sing this projection operator theorem we may now 
obtain all the necessary equations for the fasc recur¬ 
sive algorithm. Two types of recursive equations are 
of lnceresc. First, equations are needed that provide 
the m+l£t order prediction error vectors in terms of 
che mill order errors. These are called order update 
equations. Second, equations are needed chat enable 
us to update the prediction errors as a new data point 
becomes available. These equations are referred to as 
the time update recursions. 

These two sets of equations are derived below. 


V. ORDER UPDATE RECURSIONS 

In this section the order update equations for £*, 
fY, d^, d- v , u and v are derived by making use of che 
projection operator theorem. 

Consider first the error vector .f X B+l n associated 
with the optimum m+l*i. order autoregressive co¬ 
efficients. Here, m can cake on any value In the 
range 0. 1, .... p-1, where p is the desired auto- 
rectesslve coefficient order. From equations (A) and 
vilb), we see chat 


.X _ c 

t , - Prm x 

“ur+“l, n Xi —n 

(36a) 

where X-X , * (X ! S ar+ " 1 x 1 

n.aH-i n.ra • —n 

(36b) 

Y-Y , - [Y : s"* -1 ? l 

n.nn-1 n,m . Ml 

(36c) 


Applvlng the projection operator theorem to (36a) with 
A*X, B * Y, a ■ S®*' 1 x, and b-S ra+i -y , we have 

,x _c 

f m+l,n P XY - 


- P 


(S° rt ‘ L x)((S BH ' 1 y) i ?^ Y (s' ,T+l x) | \s°'* L y)''p£ Y x 


-o+l, n n 


^m,n 


m,n 


‘XY 
(37) 


As mentioned earlier, to implement the recursive al¬ 
gorithm we only need che error element at time n, that 

1® f?., (n). From equation (37) we see that 

erl,n , 


sr+l,n 


(n) 


r (n) - d* (n) 

m n y w m n 
m,n 


(38) 


Equation (38) is the desired order update equation 
for f*. 

In a similar manner, the order update equation for 
fY is found to be 


^m+l,n 


P ?X*u 


where X and Y are defined in equation (36). 
the projection operator theorem yields 


f* - f y 

-o+l,n -m, 


T m,n 

v 

m,n 


d* 


(39) 
Applying 

(40) 


The n£il component of equation (40) is the desired 
order update equation for the forward y prediction 
error, chat is 

n (n) * f m n (n) " d m n (n) (41) 

□H* i,q to • n -'in, n tn, o 


The order update equations for the delayed backward 
prediction error vectors may be similarly derived. 

These equations are, however, not as useful as the 
combined order and time update equations. The combined 
update equations give d* . . and d? . , in terns of 

—®*rl, orj. tOTx a 

x y 

d _ and d In deriving the combined order and time 

—m, n —to , n 

update equation for the delayed backward Xq estimate, 
we first note that from equation (17b) 


d x 

-m+l,n+l 

■ p h< s ^> 


(42a) 

where 


Ha » 



X • X_, 
n+1, 

m 

,aH-l • 

X 

... j 

(42b) 



L n,m. 

- n j 

i 

? ' y 

y m l 



(42c) 


n,m( -^n 


We now apply che projection operator theorem to 
equation (42a) with A“X and B«Y in equation (34). 
After some simple algebraic manipulation, we get 


-o+l,n+1 


r.o.: 

[c 



(43) 


m,n 


The n+l£t component of equation (43) yields the desired 
update equation 

u* 


Vl,n+l (n+1) 


Qk* Q 


m,o f x 
* 1 


m,n 


(n) 


(44) 


m,n 


The delayed hackward time and order update equation 
for Zq Is derived in a similar manner. The details are 


omitted, but it is readily shown that 


Cl.n+l (n+1) ' d m,n (n) 

- f y (n) 

u_ m*n 

(45) 


m,n * 


Finally, Che order update 

equations for the 

scalars 

Um.n and v m>n are derived. 

From equation (32) 


u m+l,n " £ n P X? m, 


(46a) 

where X * X . * (X I 

n,nH*l n,m . 

(S~V)1 
~n 

(46b) 

Y-Y * (Y 

n,m+l n,m . 

(8^)1 

(46c) 
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Applying the projection operator theorea glvea the 
desired result 

*Wl.n ‘ * +P XY i 


u * u 

»H t n m.n 


Also, Cron equation (33) 


Vl.Crt-l * <* 8 > 

where X and Y are defined In equation (46). Applying 
the projection operator theorea yields the coablned 
order and tine update equation for v. 


art-1, rrt-L a,n 


VI. TIME UPDATE EQUATIONS 

The remaining recursive equations update the forward 
and delayed backward errors as a new data point Is 
obtained. For this reason these equations are called 
time update equations. 

When a new daca point becomes available, the effect 
on the prediction error vectors is to append a row to 
the botton of their defining matrix equation [see, for 
example, equation (7)]. Appending a row to the 
ontrices X n n and Y n>m does not seem to fit In the 
framework of the proj’ectlon operator theorem. In which 
columns are appended to X n>m and Y n>0 . It curns out, 
however, that we can accomplish the cask of annihilat¬ 
ing a row In the error vector matrix equation by 
appending to X„ jm and Y n>m the pth basis vector defined 
by ’ ! 

^ - [0 0 ... 0 1] (50) 

To see how this works, let us consider as an example 
the forward prediction error vector for x^. If we 
append the n*l vector e„ to X nm In equation (7) we 
have 


fcn (1 >' 


f* (n) I 
l«,n J 


x(n-l) x(n-m) 


m.n < 


Equation (53b) Is satisfied only if we force 
f x (n) ■ 0. This can always be done because the 

0« u 

scalar £ appears only in che last row of equation (51). 


In particulars f£ tn (n) - 0 if we chooae 
m 

C - - I a (i)x(n-l) (54) 

1-1 

Since $jg t n(h) - 0, equation (53a) la seen to depend 
only on the first n-1 components of the vectors. It la 

easily saen, then, that f* (k) for k-1,2.n-1 are 

■»“ 

determined using the first n-1 rows of (51) In such s 

manner that the vector [f* _(l),...,f* (n-1)]' Is 

oun 1(0 

orthogonal to the first n-1 components of each column 
of Y d>b . But this is exactly the problem of determin¬ 
ing the forward prediction error baaed on n-1 data 
points. Thus, we see that 


4. -PH 


Similar arguaMnta show that this time annihilation 
property also holds for „• „» and d^ q with re¬ 

sulting formulas similar to equation (55). We finally 
note that the scalars Om,n. x B>n , u B n» •“* v m,n » r » 
formed as Inner products of the prediction error 

vectors. Since the lasc element of f* (or f y d* , 

m • o ~B,D l|D 

) i s zero, It follows that 

—m.n 


Vn-l ’ 


and similarly for x B>n , p b , and v B>n . 

With these thoughts In mind we are in a position to 
derive time update equations for a, t, u, end v. First, 
let us define the augmented matrices 

5 -‘*n,«:*nl (57 *> 

? < 57b > 

Then It follows that 

Vn-l ' £.nl + ^.n J ’ ^ P Sy (S “*V < 58 > 

where X and Y are defined by equation (57). Application 
of Che projection operator theorem yields 

°m,n-l “ inl P XY * P XY V^Xy]^ 

a m,n-l * V. - f -,n* (n)[1 ' Vn’Xn^ (59 > 


where 1 - y_ _ 4 « + P™ e 
m.n - —n XY —n 


By rewriting equation (59) we arrive at the desired 
time update equation 

( f L(n)]tdl! „WI ( 61 ) 

m.O_BjQ_ 


where f* Is used instead of £ x ro Indicate Che presence 
of the gn vector. The optimal f£ >n vector la given by 

^ n ‘ P K V (52a) 

where ^ B,a ^ 

* ' (X n.m i S„1 < S 2b) 

Y *f Y n, ra ;Vl < 52c) 

From equation (27b) we know chat 

« ' f =,n • ^ » * 0 • 1*1. 2 .» <53a) 

< fj „ , e_ > - 0 (53b) 


0 - 0 , 4 

m,n m,n-l 


The time update aquation for x Is found In a similar 
mannar by using 3 for A, ? for B, S^ljjjfor a., and ^ 
for b In equation (34) to yitld 

n< n) l*t £ « „(“>) (62) 

T . T + -StS_StS_ 

m,n m,n-l 1 - Y_ „ 

m.n 

Tht update aquation# for u and v art found to ba 


u » ii , 

m.n m,n-l 


V — V , 

m,n m,n-l 


ir (n)]*[f' (n)) 

m.n_m.n_ 


[d x (n))*[d Y (n)] 
m.n m.n 
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by using arguments analogous to those used In deriving 
equations (61) and (62). 

Finally, update equations for Y 0 ,n ate needed. To 
obcaln an order update equation, ve note that 

i — v • e * P—— e 
Y m+l,n tnXY =n 


(65a) 


where 


Y ‘ t Y n,mi 


(65b) 

(65c) 


Applying the projection operator theorem, 

1-y *e + P c e 

'rofl.n 4 

(s^x) t ^ 


1 - Y, 


m, n 


[<£ Jn)r[d> (n) ] 


m,n 


or, equivalendy 


[d* (n) 1 [d' (n)) 

m.n m.n 

y ay 4 > ■■ t— ■. —. I . 

m+l,n m,n * 

u 

m, n 


( 66 ) 


(67) 


The (undelayed) backward errors are defined by 


b 

-P,n 


-P.n 


S p x + [S°x : S 1 * 


P-1 


*%* (s ^ ; s v : 


: s 1 

: sP' 1 


4 li 


(71) 

(72) 


where a. and c. are p*l autoregressive coefficient 
vectors) By comparing equations (71) and (72) with 
equations (14) and (21), It Is readily seen that the 
(optimal) backward and delayed backward error vectors sre 
related by r 


and 


d y - 


b X , 
—m,n-l 


.9... 


a, n-1 


(73a) 


(73b) 


The n th components of equations (73) yield the 
desired results 


<r (n) 

m, n 


, (n-1) 


m,n-l 

d y (n) - b y ,(n-1) 
in n m n*l 


(74a) 

(74b) 


To obtain che time and order update equation for Y 
we note chat 


1 - v m+ 
where 

X - X 


• e 


n+1, m+1 


-n+1 

(68a) 

1 

.* i 

(68b) 

n,mj 


’ '1 
' 1 

(68c) 

n,mj 



[?..• i' 

f n+l,m+l ” ". Y 


Since the first element of e^ is zero, it follows 
chat 

1 - Y ra+l,n+l ■ ^ (69a) 


where 


X • [X - . x ] 
n,m.—n 


(69b) 

y * (Y n,mi*nl < 69c > 

Applying the projection operator theorem yields 


1-Y 


ra+1, n+1 


(f x (n)*[f y (n) ] 

1-y - n ' n a > n 

m.n * 


f t»+l,n+l Y ra,n + 


m.n 


(f„ (n)] [f Y (n) ] 

m, n to. n 


(70) 


Thus, using equation (74) we can replace the delayed 
backward error terms in che recursive equations by 
backward error terms. 

It Is also helpful to replace v by another scalar 
defined by 


o - [b x ]*[b y ] 

m,n —m.n —m,n 


(75) 


By considering Che defining equation for v, equation 
(33), along with equation (73) It is a simple mstter 
to show that 


W , * V 

m,n-l m,n 


(76) 


Table 1 summarizes the complete set of update equations 

with d x , d 5 *, and v replaced by b x , b^, and w. 

The initial conditions for the various quantities 

In Table 1 are obtained by considering the defining 

equations for these quantities. For example, from 

equation (24) we see chat If m - 0 then f? ■ x or 

-vJ,n -a 


f 0n (n) - x(n) 

From equation (71), it also follows that 
b x n ( n) - x(n) 

Similarly, we can show that 
f y n (n) » y(n) ■ x(n-q) 

b~ (n) » y(n) * x(n-q) 
u, n 


(77) 

(78) 

(79) 

(80) 


These six time update equations (61), (62), (63), (64), 
(67), and (70) along with Che six order update 
equations (38), (41), (44), (45), (47) and (49) obtain¬ 
ed in che last section comprise che fast recursive 
algorithm- 


VII. SUMMARY OF THE ALGORITHM 

At this point all necessary equations for implement¬ 
ing che fast recursive algorithm have been derived. In 
chis section initial conditions are discussed and a 
procedure for implementing the recursive equations is 
given. 

The implementation of che algorithm is conceptually 
easier If we replace the delayed backward error 
quantities by (undelayed) backward error quantities. 


Furthermore, since v n is the zero vector for n^q, 
it follows from equations (30) - (33) that 


O.n 


t. • u„ * w- “0 for n<q 
u,n u, q u,n 


(81) 


It is clear from Table 1 that no order update 
equations exist from s and t. Therefore, initial 
conditions are needed for dp in and Tp )0 for each p and 
some corresponding n. A little thought will convince 
one that fjj q +p (q+P) “ 0 by using an argument similar 

co the time annihilation argument presented in Section 
VI. Moreover, we can also show that d x , (q+p) • 0. 

p.q+p 

by using this type of argument. Thus, 


tn.q+m m,q+tn 


for m • 0,1,...,p 


( 82 ) 
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(T-l) 


Cl,n< n) 


f m,n (n) 

- 

b x .(n-1) o* /w* . 

n,n—1 o*n qi,q—L 

(T-l) 

f £fl.n (n) 


f m,n (n) 

- 

T _/4) n , 
a,n-i m.n a,n—l 

(T-2) 

Cl.n< n > 


b i.n-l‘ a -» 

- 

f x (n) t* /u* 
b«d a,n B,n 

(T-3) 

b Ll,n< n> 


n iVl) 

a, n-i 


f y (n) 0 /y 

Bf n n,n m,n 

(T-4) 

y nH*l,n 


li 

m,n 

- 

<3 T /ill , 

a,o m,n a,n-1 

(T-5) 



Ui . 

m,n-l 


p t /y 

m.n m.n m,n 

(T-6) 

a 

m.n 


a . 

m,n-l 

+ 

[f* „ ,(n-1)1/(1 - y „) 

m,n a,n-i m,n 

(T—7) 

T 

ra,n 


T m,n-1 

+ 

[b X _ (h) ] / (1 — Y „) 

a,n—l a,n o,q 

(T-8) 

u a, n 


u m,n-l 

+ 

[f x (n)]*[f^ (n)1/(1 - y ) 

m, n a, n a, n 

(T-9) 

jj 

m.n 


“m.n-l 

+ 

(b X (n)]*[b^ (n)]/(1 - y n ) 

a,n m,n m.n 

(T-10) 

Y sw-l,n 


Y m,n 

+ 

tb x „ , (n-1) ]/w* n . 

m,n-i a,n-i o,n-i 

CT—11) 

Y mel,nel 


Y m,n 


(f* .Wl*!f y n <n))/u* n 
m.n a,n m, n 

(T—12) 


where p is Che desired (maximum) autoregressive co¬ 
efficient order. 

Finally, when m « 0, 

Y 0,n * 0 < 83 > 

Although other initial conditions may be obtained, 
these initial conditions are the only ones needed to 
implement the algorithm. 

The Implementation of the update formulas can be 
divided Into three parts. First, for n^q the vector 
v,, is Che zero vector, so no operations are performed. 
For q+l^n^q+p+1, the maximum order m that can be 
used Is n-q -1. In this time Interval, as a new data 
point arrives not only are time updates performed, but 
also the model order is increased. For n>q+p+l, the 
model order remains at p and Che time updates only are 
performed. 

The Implementation of the algorithm for n^q+1 is 
summarized below. 

As the new data point becomes available: 

1) Set n - n+1 

2) These quantities are available from the last 
iteration: 

f* „ .(n-D.f?' , (n-l),b* ,(n-l)V „ , (n-1) 

m,n-i n,n—l ta,n—i m,n—1 


Table 1: Summary of Update Equations 

4) For each m ■ 0,1,...,min[p-l,n-q-3] find 


m,n-l’ ‘m,n-1’“m,n-l ,ui m,n-l 


3) New initial conditions: 


if n^p+q + 2. 


f 0,n (n > 
b 3,n (n) 
f b,n (n > 
b (5,n (n) 
Y 0,n * C 

'n-q-l,n-l 


x(n) 

x(n) 

x(n-q) 

x(n-q) 


T , , » 0 

n-q-l,n-l 


o ,t ,u ,u) using (T-7)-(T-10) 

01*0 m, n oi.n o,n 

f ^l,n (n) ’Cl.n (n) ’ b Ll.n (n) ’ b Ll.n (n) u,ln * 

(T-l)-(T-4) 

Vl.n u<,ing (T-U) 4 

5) For m - mln[p,n-q-2j find 

o ,t ,u ,u using (T-7)-(T-10) 

m,n m.n m.n n,n 9 

6) If n<p+q+l we need to add a filter order. 
Set m - n-q-2. 

Find: 

f ^l,n (n) ’Cl,n (n) ’Cl,n (n) ’ b L.l,n (n) u,ln * 

C T—1) — (T—4) 

Vl.n u8in ® (T - X) 

Vl.n’Vl.n u,in * (T * 5) ’ (T ' 6) 


for 

Set m 

a"0,l 

Find: 

• • • * 

0 

mln(p, 

m,n 

n-q-2) 

At th: 


0 , » T , "0 

m,n-l a,n-1 


to arrive. 

It is clear from Che above summary that 0(p) multi¬ 
plications and additions are required to update cha 
prediction errors. More specifically, in Che time up- 
dace mode [l.e., when n>p*q+2 so no filter orders 
need to be added! 14p multiplications and lOp additions 
are performed per update. In the time and order update 
mode, [i.e., when n<ptq+2|, 17p multiplications and 
13p additions are performed. It should be noted that 
this computational requirement may be significantly 
reduced by using a normalized lattice form similar to 
chat in [2j. 
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VIII. THE LATTICE STRUCTURE 


The faac recursive algorithm lands lcsalf co a digi¬ 
tal filter structure known as the lattice filter. This 
structure can be seen by considering equations (T-l) - 
(T-4). (t Is clearly seen that the filter that pro¬ 
duces the four outputs f x (u), f y (n), b x (n) and 


a.n n,n 

b y (n) is given In Figure 2, where 
x e.n 


y m, n 

a 7 • -*— 

“• n »n. „ 



Figure 2: 

Filter Realization of Error Order Update Equations 

Thus, the entire pi* 1 order filter la given by p stages 
of filter shown In Figure 2. This is depicted In 
Figure 3 



p stages 


Figure 3: 

Realization of che pi* 1 Order Filter 

In the startup phase of the algorltha [I.e., when 
the filter order Is less than che desired order) the 


filter of Figure 3 begins with one stage, and success¬ 
ive stages are added aa new data points arrive. Alter¬ 
nately, all p stages aay be in place at the beginning; 
however, tha four aulclpller sections of each unused 
stage is set to zero until that stage Is to be used. 

The lattice filter structure at Figure 3 nicely 
depicts the relationship between the autoregressive 
coefficients and che prediction errors. To see this 
relationship, let us denote che transfer functions froe 

x(n) to f x (a) and froa x(n) to b x (n) by F x (z) and 
■ t Q n,a ■ 

B x (z), respectively. Let us also define the auto¬ 
regressive coefficient transfer functions 

A(z) - 1 + a (l)!* 1 +...+a (m)z“* (88) 

B B ■ 

jl (z) - a (0) + a (l)z' 1 +...+S (s-l)j"* t ' 1 +t"* (89) 
BOB B 

where the a, and a^ coefficients are defined by 
equations (5) and (71) respectively. A little thought 
will convince one that the transfer functions are 
related by 

F*(z) - A (z) (90) 

B B 


B x (z) - A (z) 

B B 


(91) 


With this thought In mind we are now able co derive 
equations chat relate the error elements to the auto¬ 
regressive coefficients. From Figure 2 it Is clear 
that 

f* - Cl „ (n > < 92 > 

n,a B-l,n m-1, n 

b* „(“) „(») + „ . (n-1) (93) 

b ( q b— l,n b-1,0 b— 

By caking z transforms of equations (92) and (93) and 
using equations (90) and (91) we get 

V z) - Vl (z> +rl Cl,nVl (,) (94) 

*m (z) * *" 1 Cl,nVl ( * ) + Vl (l) (95) 

Furthermore, we can see from Figure 3 chat 

Aq(z) - Aq(z) - 1 (96) 

and from equations (94) and (95) that 

3 m (B) ’ 8 m,n (97) 

a(0) “ „ (98) 

B ra,n 

Given che error elements f* (n) and b x (n) for 

m,n n,n 

m*l,2,...,p we can use equations (92) and (93) to get 
oj.n and 6* n elements, then use equations (94)-(95) 
to get the a„(IO and a^k) elements for m»l,2,...,p 
and k • 0,1,2,...,m. Thus, we are able to obtain the 
autoregressive coefficients for all model orders from 
1 to p. 

Similarly, given the a p (k) and a-(k) coefficients for 
k-0,l,...,p, we can use (94), (95), (97), and (98) to 
get che a p, n and ® p ,n 4118 then the P * l sC order auto- 
regresslve’coefficlencs by working down froa order p 
co order 0. Then, using (92) and (93) we can obtain 
che desired error elements. 

Thus, using aquations (92)-(98) we are able co con¬ 
vert back and forth between the lattice error elements 
and the optimal autoregressive coefficients. 

As a final note, when q*0 (i.e., when an auto¬ 
regressive model is chosen) ic can be seen In Figure 3 
chat che cop half and bottom half of the filter are 
equivalent. In this case, only one half is needed, and 
the filter in Figure 3 degenerates to the AR lattice 
filter described in [1]. 
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IX. CONCLUSIONS 

In chis paper we have presented a recursive algorithm 
for obtaining che autoregressive coefficients of an 
ARMA model. The recursive algorithm is based on the 
prewindowed version of the high performance method of 
ARMA spectral estimation as described in Part 1. The 
recursive algorithm is computationally fast, requiring 
O(p) additions and multiplications to update tha para¬ 
meters. Moreover, the algorithm can be Implemented 
using a lattice filter structure offering numerical 
robustness and nice convergence properties associated 
with lattice type algorithms. 

We have not yet discussed the problem of recursively 
estimating the moving average coefficients in the ARMA 
model. We do not at this time have such an algorithm. 
However, it is worth noting that the moving average 
information is present in the output prediction error 
sequences, and the utilization of this information for 
moving average coefficient estimation is currently 
under study. Another area under study is the use of 
various normalization procedures Co effect a decrease 
in computational requirements and in sensitivity. 

Finally, we noce that the recursive algorithm pre¬ 
sented here is based on approximating a set of p Yule- 
Walker equations. It has been shown In Part 1 of this 
paper thac for short data lengths. Improved spectral 
estimates result from using more than p Yule-Walker 
equations. For chose cases in which the amount of data 
is small, a fast recursive algorithm based on the 
approximation of t >p Yule-Walker equations would often 
prove useful. Such an algorithm Is currently being 
pursued. 
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